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This study belongs to a series devoted to using the Szekeres inhomogeneous models in order to 
develop a theoretical framework where cosmological observations can be investigated with a wider 
range of possible interpretations. While our previous work addressed the question of cosmological 
distances versus redshift in these models, the current study is a start at looking into the growth 
rate of large-scale structure. The Szekeres models are exact solutions to Einstein's equations that 
were originally derived with no symmetries. We use here a formulation of the Szekeres models 
that is due to Goode and Wainwright, who considered the models as exact perturbations of a 
Friedmann-Lemaitre-Robertson- Walker (FLRW) background. Using the Raychaudhuri equation we 
write, for the two classes of the models, exact growth equations in terms of the under/overdensity 
and measurable cosmological parameters. The new equations in the overdensity split into two 
informative parts. The first part, while exact, is identical to the growth equation in the usual 
linearly perturbed FLRW models, while the second part constitutes exact non-linear perturbations. 
We integrate numerically the full exact growth rate equations for the flat and curved cases. We find 
that for the matter-dominated cosmic era, the Szekeres growth rate is up to a factor of three to 
five stronger than the usual linearly perturbed FLRW cases, reflecting the effect of exact Szekeres 
non-linear perturbations. We also find that the Szekeres growth rate with an Einstein-de Sitter 
background is stronger than that of the well-known non-linear spherical collapse model, and the 
difference between the two increases with time. This highlights the distinction when we use general 
inhomogeneous models where shear and a tidal gravitational field are present and contribute to the 
gravitational clustering. Additionally, it is worth observing that the enhancement of the growth 
found in the Szekeres models during the matter-dominated era could suggest a substitute to the 
argument that dark matter is needed when using FLRW models to explain the enhanced growth 
and resulting large-scale structures that we observe today. 

PACS numbers: 98.80.Es,98.80.-k,95.30.Sf 



I. INTRODUCTION 



Our modern era of cosmology has led to not only a wealth of astronomical observations, but also to two major 
conundrums, namely, dark matter and dark energy. It is therefore essential to explore theoretical frameworks that 
allow for a wider range of possible interpretations of the cosmological data. 

Such a framework can be provided by inhomogeneous cosmological models that are exact solutions to Einstein's 
equations. A good deal of theoretical work has been done on such models in the exact theory of general relativity, 
but little work has been done to compare them to observations. This is not surprising, since it is not straightforward 
in these models to derive observable functions ready to be compared to cosmological data. See, for example, fH-fl8j 
and references therein. 

This study is part of a series where we consider the Szekeres inhomogeneous models in order to develop a framework 
in which to analyze current and future cosmological observations. The models were originally derived by Szekeres 
in (l9l [20| as an exact solution with a general metric that has no symmetries and has an irrotational dust source. 
The models are regarded as the best exact solution candidates to represent the true lumpy universe we live in. For 
example, the models are put in the same classification as the observed lumpy universe in 12 ill . The models have been 
investigated analytically and numerically by several authors; see, for example, 0, [IH, |22T - |36| . 

Any cosmological model must pass at least two types of cosmological tests. The first is related to measurements 
of the expansion history and cosmological distances, and the second is related to measurements of the growth rate of 
structure formation. In previous studies fl6l |35|. we explored distances versus redshift in the Szekeres models, while 
the current work looks into the growth of structure. 

In this paper, we study the growth rate of structure in flat and curved cases of Class-I and Class-II Szekeres models 
using a reformulation of the models that was introduced by Goode and Wainwright [25|, [26[ , where the models can 
be considered as non-linear exact perturbations of the Friedmann-Lemaitre-Robertson- Walker (FLRW) background. 



* Electronic address: mishak@utdallas.edu 



2 



This formulation is well-suited for the study of the growth and structure formation in these models. The work of 
Ref. [39| extended the flat case of Class-II in J25| to include a cosmological constant and a discussion of the growth. 
Our work here extends that of [25[ to spatially curved cases in Class-I and Class-II in analyzing the growth. We also 
express the growth differential equations as well as equations for the shear and tidal gravitational field, all in terms of 
the under/overdensity and measurable cosmological parameters, which offer further insights over the metric functions 
themselves. We analyze the time evolution of the growth factor, shear, and tidal gravitational field scalars in fiat 
and curved cases and discuss their interrelationship via the propagation equations to produce stronger gravitational 
clustering in Szekcrcs models. 

The paper is organized as follows. After presenting the formalism in section II, we show in section III how the 
Raychaudhuri propagation equation gives exact nonlinear growth equations in terms of the under/overdensity that 
divide into two meaningful parts. We express these equations for the flat and curved cases of Class-II in terms of 
growth factor, the scale factor, and measurable cosmological parameters. We perform numerical integrations of these 
new equations and plot the results for various cases. The results are also compared to those of the linearly perturbed 
Einstein-de Sitter model and the nonlinear spherical collapse. In section IV we repeat the analysis for models of 
Class-I. We then explore the shear and the tidal gravitational field and their relation to the growth in section V. We 
conclude in section VI. Throughout the paper we use units such that 8nG = c = 1. 



II. THE SZEKERES MODELS IN THE GOODE AND WAINWRIGHT REPRESENTATION 

Since the original derivation of the models by Szekeres in [2(| , they have been reformulated in at least two other 
different sets of coordinates. The second formulation was proposed by Goode and Wain wright [ 25l |26| . and a third set 
was used in e.g. [l,|36j]. The relationships between these formulations can be found in [^.l26ll36|. As mentioned earlier, 
we chose in this study to use the representation of Goode and Wainwright because it is well-suited for the study of 
the growth of structure and exact perturbations of a smooth FLRW background. In order to be self-contained, we 
repeat and summarize here the presentation of the models in the Goode and Wainwright presentation (some typos 
have been fixed from their original paper). This also serves to set the notation to be used in the paper. The Szekeres 
metric in the Goode and Wainwright set of coordinates {t, x, y, z} reads 

ds 2 = -dt 2 + S 2 [e 2v (dx 2 + dy 2 ) + H 2 W 2 dz 2 ] , (1) 

where 

H(t,x,y,z) = A(x,y,z) - F(t,z) 

= A(x,y,z)-(J3 + f++l3-f-), (2) 

and H(t,x,y, z), S(t,z), and W{z) are all positive functions. /3 + and /3_ are functions of z only, and /+ and /_ are 
functions of t and z. A(x, y, z) and v{x, y, z) will be specified further below. 

The source of the spacetime is irrotational dust, and the coordinates are comoving and synchronous. The cosmic 
dust fluid thus has the 4- velocity vector u a = [1, 0, 0, 0]. 

The two classes of solutions are listed below with further specifications of the functions for each class. The function 
S(t, z) satisfies the generalized Friedmann equation 

*> = -* + f, (3) 

where ' = d/dt, M = M(z), and k = 0,±1. Next, it can be verified that the functions /+ and /_ (see e quat ions © 
and ([9]) below) arc the increasing and decreasing solutions of the following ordinary differential equation |26| 

F + 2p-^F = 0, (4) 

which can be derived from the field equations as well as from the Raychaudhuri equation as we discuss further in the 
following section. 

Interestingly, the time evolution of the two spatial classes of the models is fully described by the same two equations 
([3]) and (j4|). Equation ([3]) reduces to the usual Friedmann equation for a fixed z while equation (j4|) governs density 
fluctuations. 

The matter density in the models is given by 

. 6M / F\ 6MA 
p{ t, x , y ,z) = — I ! + -)=—. (5) 
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In this formulation, the vanishing of the functions /?+ and /?_ is the necessary and sufficient condition for the models 
to reduce to FLRW models. In this limit, the sign of the matter density is determined by that of M(z), and thus we 
shall limit our interest to models with M(z) > 0. 



A. Time dependence 

The time dependence of the models is specified by the parametric and implicit solutions for the functions S(t,z), 
f+(t, z), and f-(t, z) from equations ([3]) and (|4]). The solution of the generalized Fricdmann equation ([3]) is given in 
parametric form using n as: 

s = M ^Ml with t -T(z) = MhM, (6) 
dn 

where 

!n — sin 77, k = +1 
sinh77-?7, k = -1 (7) 
n 3 /6, k = 0. 

(Note that M(z) > is required in all cases, and that when fc = or = -1 then we require that S > 0.) So the 
scale function S has the same time dependence as that of an FLRW dust model. Next, the solution to equation (|4]) 
gives 

(6M/S) [1 - (77/2) cot(r//2)] - 1, k = +1 
f+={ (6M/£9[l-fa/2)coth(»7/2)] + l, k = -1 (8) 
r/ 2 /10, k = 

(6M/S) cot(?7/2), k = +1 
/_ = { (6M/S) coth(r ? /2), k = -1 (9) 
24/r/ 3 , fc = 0. 

B. Spatial dependence 

The spatial dependence of the models separates them into two classes. 

Class-I: S = S(t, z), S z ^ 0, /± = f±(t, z), T = T{z), M = M(z), 

with 

e v = f(z)[a(z){x 2 + y 2 ) + 2b(z)x + 2c(z)y + d(z)}- 1 , (10) 

ad - b 2 - c 2 = e/4, e = 0,±1, 
A = fv g -k0+, W 2 = {e-kf)-\ 
/3+ = -kfM z /(3M), /3_ = /T,/(6M), (11) 

where the subscript z in the equations means differentiation with respect to z. 

Class-II: S* = 5(f) , f± = f±(t), T = const, M — const, 

with 



e 



v = [l + k/4(x 2 +y 2 )}-\ k = 0,±l, W = l, (12) 



A _ , e»{a(z)[l- ±(x 2 +y 2 )]+b(z)x + c(z)y}~k(3 + , k = ±1 r * 

1 "( 2 ) + 6(z)a; + c(z)y-/3 + (x 2 +?/ 2 )/2, fc = 



(7 



4 



C. Raychaudhuri equation and gravitational attraction 

The complete dynamics of a cosmological model can be expressed in terms of a set of evolution and propagation 
equations (see, for example, (2ll |38| and citations therein). One of the fundamental propagation equations is the 
Raychaudhuri equation [37j , and it can be viewed as the basic equation of gravitational attraction pTj . In the case of 
the Szekeres irrotational dust, the equation reads 

e + ie 2 + 2<r 2 + ^ = o, (14) 

which includes the following quantities associated with the 4- velocity vector u a : 
• the rate of volume expansion scalar 

© = U°;„ = 3| - £ (15) 



the rate of shear tensor 



S H 



°ab = U(o;6) + "(o«6) ~ (flab + "aMt), (16) 
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where for our models the corresponding non-zero components are 

2o* x = 2a« v = -o' z = -~ (17) 

• and the matter density, p, which is given by equation ([5]). 

In the above, a semicolon denotes covariant differentiation, and parentheses around indices indicate symmetrization. 

For a complete Raychaudhuri equation including pressure, vorticity, 4-acceleration, and a cosmological constant (all 
zero in our models here) , see for example [2lL l38j . 

Now, using equations (|15p . (|17[) , the generalized Friedmann equation ([3]) (making the appropriate substitution for 
S), and the density equation ((SJ), we can write the Raychaudhuri equation (|14[) as second-order ordinary differential 
equation that is linear in the function F: 

F + 2p-^F = 0. (18) 

As pointed out in [26j j , this differential equation of the metric function F derives from the field equations as well as 
the Raychaudhuri equation. However, as we will explore in the rest of the paper, associating this equation with the 
meaning of the Raychaudhuri equation and that of other evolution equations involving shear and tidal gravitational 
fields will help us building a consistent discussion of gravitational clustering in the Szekeres models. 

Including the previous section, all equations up to this point apply generally to both Class-I and Class-II of the 
Szekeres models. However in the following sections, we need to treat the discussion of the growth equations for Class-I 
and Class-II separately in view of their different spatial dependencies. 



III. STRUCTURE GROWTH EXACT EQUATIONS IN SZEKERES CLASS-II MODELS 

We start by exploring Class-II first because of its simpler spatial dependence. In this case, the function S(t) is only 
a function of t while M and T are constants. The connection to the growth of large scale structure is simpler since p 
in equations (fT9|) and (|20[) below is only a function of t. 

Now we define, in the usual way, the density contrast S as 

6bx,v,z)= rt> X > V :$-M , (19) 

so that we can write 



p(t,x,y,z) = p(t)[l + S(t,x,y,z) 



(20) 
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Comparing equation ([20)) and equation ([5]) specialized to Class-II, that is 



6M 

Pit) = ^ (21) 

which is a function of t only since M is constant, we identify immediately 

5{t,x,y,x) = — . (22) 

Using the relations F = SH + SH and F = SH + 2SH + SH obtained from equation (|22p , the Raychaudhuri equation 
(fT5| becomes 

... S ■ M 

SH + 2SH + SH + 2-{SH + SH) - 2,—SH = 0. (23) 
b b 6 

Noting that H = — F and H = —F, we can rewrite the above equation (divided by H) as 

' s+2 §*- 3 P 2+ 4*- 3 P= o - (24) 

We can eliminate H/H with the above relation for F and arrive at the following exact growth equation for the Szekeres 
Class-II models: 

5, 9 ^(*)i o M * 2 i2 o M A 2 n ^ 

5 + 2 5(I) 5 " 3 sW ~T+s 3 sW = °- (25) 

Now, we use the same idea of Goode and Wainwright [l^, |2f|, who stated that the Szekeres models in their 
formulation allows one to compare them with linear perturbations of the FLRW models [25| . Additionally, as explained 
in (26[, the function S(t,z) in the generalized Friedmann equation ([3]) corresponds to the scale factor in the FLRW 
models in the sense that for each value of z it satisfies the Friedmann equation. In other words, for each Szekeres 
model explored, we associate a corresponding linearly perturbed FLRW model. The exact growth equation obtained 
from the Szekeres model is then compared to the linearly perturbed associated FLRW model. The unperturbed 
FLRW model is what is called here and in other papers the associated FLRW background. It is worth clarifying that 
this association is not based on the usual procedure where limits are imposed on some coordinates or special values 
imposed on the metric functions such that the Szekeres models will reduce to an FLRW model Q , even if it may be 
related to it. 

So following [ID, [2(| , we associate the Szekeres models to non-linear exact perturbations of an associated FLRW 
background, and we write here the growth equation (|25|) with an exact density contrast as perturbations of an FLRW 
model with the corresponding Hubble function and matter density as given below in equations (|2T|) and (|28[). The 
equation then reads 

5 + 2M(t) FLmv _j - 47rGp(t) FLRW _ B 5 - - 4irGp(t) FLRW _ B 5 2 = 0, (26) 

where we have identified the Hubble expansion rate of the FLRW background (FLRW-B): 

m) FLRW - B = j^- (27) 

Note that this H(t) is not the same function as the metric function H(t,x,y, z). Since M is a constant in Class-II 
here, we have used equation ([5} to express the factor in terms of the smooth background matter density as 

=^Gp{t) PLRW _ B . (28) 

We note that while throughout this treatment we have set k = 8irG = 1, we restored its value here just to make the 
identification with the usual FLRW expressions immediate. 

It is also worth noting that the Raychaudhuri propagation equation provided us with an exact equation for the 
growth function S that can be remarkably split into two meaningful parts. One part, which consists of the first three 
terms, is identical to the usual equation of the growth for a linearly perturbed FLRW model. The other part contains 
the remaining nonlinear terms in S and S and is similar to second order perturbation terms. But we note here that 
our equation is exact, and 6 does not have to be small. In the next sub-sections, we specialize the exact growth rate 
equation (f2"5)) of the Szekeres Class-II to the flat and curved cases and then integrate them numerically. 
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A. Growth rate equations and integration in flat Szekeres Class-II models 



In this case, the corresponding FLRW background is thus the well-known Einstein-de Sitter model that is appropriate 
for describing the matter-dominated phase of cosmic evolution. For an Einstcin-dc Sitter background, the scale factor 
can be solved from equation as 

/Q \ 1/3 

S(t)=(-M) t 2 / 3 = i 2 / 3 , (29) 

where for k = we can set M = 2/9 in all generality (for example, see [26|). Using equation (f2"9")l in equation (|2l))) or 
(|25| yields 

' + ^-^ s -ih" 2 -P 2=0 - (30) 

We note again that the first 3 terms are exactly the same as the ones in the usual linearly perturbed Einstein-de 
Sitter case, so this could be identified with the linear theory when <5 is small. The second part is made of two nonlinear 
terms in 5 and 5, and can be compared to second order perturbation terms [40j . 

The full equation obtained can be compared as well to the spherical nonlinear collapse model, and we do that further 
below. As in the case of nonlinear spherical collapse (40| for Einstein-de Sitter, a parametric solution is well-known 
(see, for example, (ij), but we found it more practical for our cosmology numerical codes to perform the integrations 
numerically. Moreover, the numerical integration schemes are expandable to cases where parametric solutions are not 
known. 

Next, we write the growth rate equation (|30p in terms of the scale factor. To convert time derivatives to S derivatives, 
we first note that 



and 



i = 3^', (3D 

l = vs s '-w s "- < 32 » 

where primes denote derivatives with respect to 5*. Here and throughout the rest of the paper, we write a(t) for the 
scale factor S(t) by analogy with FLRW models and take it to have the standard normalization a(to) — 1 today. We 
thus obtain 

2 a 2 a 2 1 + 5 2 a 2 vy 

For our numerical integrations, we write the growth equation in terms of the growth rate G = 5 /a. Using the relations 
5' = G + aG' and 5" = aG" + 2G', we get 

7G' 2 {G + aG') 2 3 G 2 
G + 2~~a 1 + aG ^7"°' (34) 

We integrate this equation using a fourth-order Runge-Kutta algorithm with adaptive step size [52] ■ We implement 
the Runge-Kutta code with the function vectors y = {G, G'} and 4^ = {G', G"} [42[ so that the second-order ODE 
(|34|) is reduced to two first-order ODEs that are immediately integrable. 

Our results for the integration are given in the left part of figure [T] The horizontal line G = 1 corresponds to the 
linearly perturbed Einstein-de Sitter model, while the red curve (highest) represents the solution to the full exact 
growth equation (j34|) for the flat Szekeres model of Class-IL The Szekeres growth is up to a factor of 3 stronger than 
that of the perturbed Einstcin-dc Sitter background. 

For further comparisons, we also integrate numerically the growth for the well-known nonlinear spherical collapse 
model in an Einstein-de Sitter background (see, for example, [4(|), where the governing equation is given by 

A"^ 3<5 ' 36 4 1 A' 2 3 62 n 



or in the G-notation 



7G' 4 (G + aG'f 3 G 2 
G + 2~'Ya 1 + aG ~ 2~ °" (36) 
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The resulting integration is also plotted in the left part of figure 1. The Szekcres growth is found to be stronger than 
that of the spherical collapse model, as well. Moreover, the relative difference also increases as a function of time. 
This indicates that the growth rate is different when the spherical symmetry approximation for inhomogeneities is 
not assumed. As we discuss in section V further below and in the concluding remarks, the Szekeres models have shear 
and tidal gravitational fields that contribute to the enhancement of its gravitational collapse and growth rate. 



B. Growth rate equations and integration in curved Szekeres Class-II models 

We derive now the growth ODE for the positively and negatively curved cases (k = ±1). Both situations can be 
handled simultaneously and lead ultimately to the same equation that we integrate numerically. We begin from the 
growth equation (|25|) 

6 + 2 *5-3^5-^-6 2 -3%5 2 = 0, (37) 
a a 3 1 + 5 a 3 y ' 

and the generalized Fricdmann equation ([3]) 



2 



2M k 



(38) 
a 6 er 



where an overdot denotes differentiation with respect to t (partial in the case of 5 and total for a). We define 
Qm = 2M/(a 3 H 2 ) and flk = — fc/(a 2 H 2 ) by analogy with standard cosmology such that (|38|) can be written as 
SIm + Slfc = 1. We proceed with a more general approach than for the flat case that does not require an explicit 
expression for a as a function of t for the different k values. This is accomplished by using equation Q38p directly to 
recast equation ([37]) into integrable form. As before, we convert time derivatives to a derivatives via the (now general) 
relations 

• 85 85 da ,, 

d = Yt = 8-aTt =5a (39) 

and 

5 = %-(6'a) = a 2 5" + M'. (40) 
at 

Substituting these relations into equation (|3"T|) . we find 

M „ 2 „, „ M 



5" + ^ + - )5' - 3^5 - -^-5' 2 - 3^5 2 = 0. (41) 
\a z a J a a a z 1 + 5 a^a^ 

We next use equation Q38p to eliminate a and a in favor of Qm ^ n d (implicitly) £1/.. Taking its time derivative and 
dividing by a 2 gives the coefficient on 5': 



a 2 4-fi 



M 



(42) 



(43) 



a 2 a 2a 
The coefficients on <5 and 5 2 are also easily converted as 

3 m = zn_ M 

a 3 a 2 2 a 
and so we obtain the equation 

„„ / n M (a)\ 5' 3 , s 5 2 , 9 3 , J 2 

5" + (2 - -f±) - - -n M[a) - - - -n M (a)- = 0. (44) 

Again, we see that the first part agrees perfectly with the usual growth equation for the spatially curved linearly 
perturbed FLRW models , while the second part is made of two nonlinear terms that can be compared to second 
order perturbation terms. Converting as before to G-notation for integration, we obtain 

VM(a)\ G> , n/1 n G 2 (G + aG') 2 3 G 2 

a z a 1 + aG 2 a 



G" + ( 4 - ) - + 2(1 - n M (a))^ - ",r: r< > ^ M (a)— = 0. (45) 
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FIG. 1: LEFT: Growth rate of structure in flat Szekeres Class-II models (or Class-I models with a fixed value of z; see section 
IV) (solid-red curve), the usual spherical collapse model (green-dashed), and the perturbed Einstein-de Sitter (EdS) model 
(blue-dotted). The Szekeres growth rate is stronger than that of the perturbed EdS by up to a factor of 3. The Szekeres 
growth rate is also stronger than that of the spherical collapse model. CENTER: Growth rate of structure in positively curved 
Szekeres Class-II models for various values of Q]^ (or Class-I models with a fixed value of z and various values of Sl° M {z)). The 
growth rate in linearly perturbed FLRW models with the same values of are plotted for comparison. RIGHT: Growth rate 
of structure in negatively curved Szekeres Class-II models for various values of Sl° M (or Class-I models with a fixed value of z 
and various values of €l° M {z)). The growth rate in linearly perturbed FLRW models with the same values of Q^j are plotted 
for comparison as well. In both cases, the Szekeres growth rates are stronger than those of the corresponding perturbed FLRW 
models by up to 5 times. 



We must now express J7m explicitly in terms of the scale factor, a, and the matter density parameter today, 
Using the definitions of D,m and Q,j- and the Friedmann equation (|38|) evaluated today, we see that 

1-^ = ^ (46) 
allows the Friedmann equation (|38p at a given time to be written as 

H = — - ?±J£ L 47 



T2 



where a super- or subscript naught means that the parameter is evaluated today. Dividing by H , (|47|) becomes 



and dividing (J3HJ) by H 2 , yields 



Thus we find the usual relation 



l-U M {a) , (48) 



t = ±[SP M -a(Sl° M -l)]. (49) 
u a 



""M- m + ff-q,) - (50) 



G" + |4-..„„ ^g +2 fl-^ «*■ „„,U-? (G + °?'-fL. ^g=0. 



which, upon substituting into P5|) . gives finally 

nSf \9L +2 ( 1 

2[Q.° M + a(l-il° M )]J a \ £l° M + a(l - Sl° M ) J a 2 a 1 + aG 2 \n° M + a{l - fl° M ) j a 

(51) 

It is this equation that we integrate numerically in our code for various values of Q° M - We note that the same 
differential equation (|51l) governs both the k = +1 and the k = — 1 cases. The dependence on f2fc is accounted for in 
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FIG. 2: LEFT: Energy density in flat and positively curved Szekeres Class-II models (or Class-I models with a fixed value 
of z; see section IV). The density is plotted as a function of the scale factor well within the matter-dominated cosmological 
era (a < 0.35) (i.e., well before a w 0.60, the scale factor associated with matter and cosmological constant equality in the 
FLRW-LCDM standard model [13, (3) • The density diverges as it should as a tends to the initial singularity, and it decreases 
monotonically for increasing a with no further divergences during the matter-dominated era. CENTER: Energy density in 
negatively curved Szekeres Class-II models (or Class-I models with a fixed value of z). The flat case is also given as reference. 
The same observations hold here. Further discussion is given in sections III-B and V. RIGHT: Comparisons of some Szekeres 
models to a growth range band of 5% around the Einstein-de-Sitter growth, which is also the approximate representation for 
the FLRW-LCDM model during the matter dominated era under consideration. We find that there are Szekeres models that 
are consistent with such a growth range during the cosmological era considered, but the Szekeres models require only about 
a third of the matter density compared to that of Einstein-de Sitter. This is consistent with our result of stronger structure 
growth in the Szekeres models. 



the sign of 1 — f2°y. As expected, it is clear that setting Q° M = 1 indeed recovers the growth equation for the k = 
case derived in the above sub-section. 

Our results for the integration of the curved cases are given in figure 2 for various values of ft^.j and include positively 
and negatively curved models. We also integrate and plot the corresponding linearly perturbed FLRW models. As 
in the spatially flat cases, we find that the curved Szekeres growth is up to 5 times stronger than that of the linearly 
perturbed curved FLRW. Also, some Szekeres models with f2°/ < 1 can still have a stronger growth than that of the 
Einstcin-de Sitter growth. 

In addition to the growth, it is worth plotting the energy density evolution in order to verify that, after the expected 
initial singularity, it has no other divergences during the matter-dominated era of interest here. We use equations ([5]) 
and (|20[) and some of the steps above to write the energy density as follows: 

p=^±(l + 6)=3M 2 a -Mf(l + aG). (52) 

Our plots in figure [2] for the density as a function of the scale factor for the flat and curved Szekeres show that the 
energy density diverges toward the initial singularity, as it should, and after it decreases monotonically with no other 
divergences during the matter-dominated era plotted here with a < 0.35 (i.e., well before a ~ 0.60, which corresponds 
to the equality time between matter dominance and cosmological constant dominance in an FLRW-LCDM (Lambda- 
Cold-Dark-Matter) model; see, for example, [43|,[44|). Our first goal here is to represent the growth rate during the 
matter dominated era with a < 0.35, and as the density plots show, we are far away from any possible pancake 
singularity that could occur at later times (2(| . Additionally, such pancake singularities can be avoided in some cases 
by the addition of a cosmological constant to the models [39||, although the discussion there was limited to the flat 
case. Finally, we conclude in this section that the growth rate of large-scale structure is consistently found to be much 
stronger in the Szekeres models than in the linearly perturbed FLRW models with the same matter density. 
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IV. STRUCTURE GROWTH EXACT EQUATIONS IN SZEKERES CLASS-I MODELS 



A. Flat Szekeres Class-I models 



Unlike in the spatially curved Class-I models, the time evolution equations ((6]) in the flat case can be decoupled 
and one can set without loss of generality M = 2/9 (26[. As a result, the treatment of the growth in the flat case can 
be developed in exactly the same way as was done for Class-II in section III-A. However, a closer look at equation 
(fTTT) shows that for k = in Class-I, (3 + = 0, so there are no growing modes and this case becomes of no interest for 
our purpose. 



B. Growth rate equations and integration in curved Szekeres Class-I models 



In the spatially curved Class-I models, the functions S(t,z) and M(z) are now z-dependent. As a result, the 
identification of an under- or overdensity function as well as a background will require some discussions and definitions 
different from the ones we made for Class-II. 

We define a density contrast 5 as 

xn \ _ P(t,n,y,z) - p q (t,z) 

o{t, x, y, z) = , (53 

Pq{t,Z) 



where we use the quasi- local average density [45H48I ] : 

x SzSxJ y J r p(t,x,y,z)V=hdzdxdy 
P ^ Z) = UJ y F^hdz dxd y = W^M* (54) 

for the domain T>[z] delimited by z ^constant and contained in the hypersurface £=constant. Here h is the determinant 
of the 3-dimensional part of the projection tensor h a t, = g a t, + u a ui,, and J 7 is a physically meaningful weighting function 
that was discussed in I45rj49l |. 

It was shown in that p q is a coordinate independent quantity that can be expressed in terms of curvature 

invariants and is given for the Szekeres Class-I models by [48| 

6M(z) 

p ^ z) = w^- (55) 

As discussed in [45l - |47j , the term quasi- local follows from the integral definition used for the quasi- local Misner-Sharp 
mass-energy in spherical symmetry (50j . Bearing in mind that M(z) is associated with such a quasi-local mass-energy, 
the result (f5"5j) is not unexpected. 

Using the definition (|54[) we can write 

p(t,x,y,z) = p g (t,z)[l + S(t,x,y,z)]. (56) 
Comparing equation ([55]) and equation ([5]), we identify 

F 

5(t,x,y,z) = — . (57) 

Again, using equation ((57)) in the Raychaudhuri equation (|T8)l gives, after a few steps, the following exact growth 
equation for the Szekeres Class-I models: 

$ , J(t,z) i M(z) . 2 ~ 2 M(z) S2 

6 + 2 w^) s " 3 W^r d ~ m " 3 W^r 6 = °- (58) 

Next, following a large number of papers using the Lemaitre-Tolman-Bondi models (LTB) for the purpose of comparing 
them to observations (see for example @, BUll, E^j and references therein), we can generalize some definitions from 
the FLRW models to the Szekeres models. We start by rewriting the Szekeres generalized Friedmann equation ([3]) as 



(59) 
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We then define 

0,,(f r\ = 

S(t,z) 3 H(t,z 

and 



n u (t,z)= ou 2M S z l „ (60) 



n " (t ' g)s 5(MWM) a (61) 
by analogy with standard cosmology such that (|50[) can be written as 

n M (t,z) + n k (t,z) = i. (62) 

Now, following the same approach we used for Class-II earlier and the idea introduced by Goode and Wainwright 
pM |26| . we can define a background that is independent of coordinates x and y and where the S(t, x, y, z) can be 
regarded as an exact perturbation satisfying the exact equation 

<5 + 2H(t, z)6 - 4irG Pq 5 2 —l - AirG Pq S 2 = 0, (63) 

1 + 5 

Again, we can observe that this exact growth equation splits into two meaningful parts. One part, which consists 
of the first three terms, is similar to the usual equation of the growth for a linearly perturbed FLRW model. The 

other part consists of the nonlinear terms in S and S. 

We follow now the steps performed in the curved cases for Class-II with an overdot here always denoting partial 
differentiation with respect to t. We use equation (|5i?|) to recast equation ((55)) into intcgrable form. After some steps 
similar to equations (|39[) - (|43[) and moving to a-notation, we obtain 

2 -n M (a,z)— r(5 - -fi M (a,z)— = 0. (64 

2 J a 2 a 1 + 5 2 a z 

It is worth noting again that here a is a function of t and z while it was a function t only in Class-II. 
Next, converting as before to notation using G = 5/a for integration, we obtain 

~„ / f2jvf(a, z)\ G' n , i , ..G 2(G + aG') 2 3 , ,G 2 n 

G" + U \2-L — + 21- Q M (a, z))— - J- - -Q M (a, z)— = Q. 65 

V 2 J a a 2 - a 1 + a G 2 a 

Now, we can express fijvf(a, z) explicitly in terms of the scale factor, a, and the matter density parameter today, 
CIm{z) where a super- or subscript naught means that the parameter is evaluated today. One obtains 



!lMlM, -n»( 2 )|a(l-(!»( 2 )) : 



(66) 



Mr/ ^ u y ± M 

which we substitute into (p5|) to finally get 



Sl° M (z) \ & 2 A OS,(z) ^ G 2(G + aG') 2 3^ ^ G 2 



2[OS,(«) + a(l a V + « 2 « 1 + aG 2 + a(l - , 

(67) 

We first note that unlike the case of Class-II, here tt° M (z) has a z-depcndance and in order to perform the plots 
we will use it for a fixed value of z. Just as we specified for the quasi- local average density, we use this quasi-local 
average density number within a domain limited by such a constant value of z. 

Our results for a fixed value of z are then identical in form to those of Class-II and are plotted and discussed in 
figures 1 and 2. The comments and remarks made for Class-II apply here as well. 

In summary, for both Class-I or Class-II, the Szekeres growth is found to be up to 5 times stronger than that of 
the linearly perturbed curved FLRW. Also, our plots in figure [5] for the density as a function of the scale factor for 
all classes and the cases that we explored show that the energy density diverges toward the initial singularity, as it 
should, and afterward it decreases monotonically with no other divergences during the matter-dominated era that we 
considered in this analysis. In order to explore further the strong growth found in the Szekeres models, we examine 
in the next section the shear and gravitational tides in the models. 
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V. THE SHEAR AND THE TIDAL GRAVITATIONAL FIELD IN SZEKERES MODELS 

In this section we analyze the time evolution of further physical quantities that can affect directly or indirectly the 
gravitational clustering and the growth rate of structure in the Szckercs models. The source being irrotational dust 
(i.e. zero vorticity and 4-acceleration) , the next quantities of interest are the rate of shear and the tidal gravitational 
field. For this, we consider their corresponding scalar invariants. 

The discussion in this section applies to Class-I and Class-II since the evolution equations and scalar invariants are 
common to both classes. The distinction between the two treatments is that when we talk about Class-II, we use 
5(t, x, y, z), G = 5/a, M, a(t), and fij^(a), while we use for Class-I 5(t, x, y, z), G = 5/a, M(z), a(t, z), and ^M( a , z) 
evaluated at a fixed value of z, as defined in the two previous sections. For simplicity of notation, we use here the 
formalism as used for Class-II but with observations and conclusions applying as well to Class-I with z fixed. 

First we consider a shear scalar that is equal to twice the commonly used magnitude of the rate of shear [111, HH 
and defined as 



c b<y a, 



(68) 



where the shear tensor components are given in equation (|17|) (note there is no factor 1/2 in the definition). Computing 
this and taking the square root, we find 



o = \ - 



2H 
~3H' 



(69) 



In order to analyze the time evolution of this scalar, we express it as a function of the scale factor and current 
measurable cosmological parameters. We can relate it to the density contrast via equation (|22j) and by noting that 
due to the metric function dependencies, we have H = —F. Differentiating 5 with respect to time allows us make the 
identification 



H 
H 



1 



(70) 



Now, using equation (|39|) and converting undcr/ovcrdcnsitics (5s) to growth rates (Gs) as before, we obtain the shear 
in terms of the scale factor and the Hubble parameter H: 



2 G + aG' 



3 1 



iG 



aM. 



(71) 



We use equation (|4"9")l to rewrite the Hubble parameter H in terms of its measured value today, Ho, and that of the 
matter density today, 



2G + aG' 



3 1 



iG 



oo 

M 



(V°m ~ !)■ 



(72) 



Along with \a\, another related quantity that is worth analyzing in our models is the shear scalar over the expansion 
scalar, JgL This is particularly useful when looking at the early time divergences as we show further below. Using 
equation (|15[) and after a few steps we find that 



e = 



3 - 



(G + aG')a 



1 + aG 



oo 



3 - 



(G + aG')a 
1 + aG 



(73) 



so that 



H 

e 



a(G + aG') 
3 3 + a(2G-aG')' 



(74) 



Before plotting the time evolution of these scalars, we look again at the Raychaudhuri equation (|14[) where one can 
view the shear scalar, er 2 , as an "effective source" that contributes along with the energy density, p, to cause stronger 
gravitational attraction and collapse. The equations and ([701 for a 2 above show that indeed the shear introduces 
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FIG. 3: LEFT: Shear in positively curved Szekeres Class-II models (or Class-I models with a fixed value of z). The shear 
scalar \a\ and scalar ratio \o\/Q are plotted as functions of the scale factor well within the matter-dominated cosmological 
era (a < 0.35) (i.e., well before a « 0.60, the scale factor associated with matter and cosmological constant equality in the 
FLRW-LCDM standard model 0, HJ). The top figure shows that the shear diverges when a tends to the initial singularity, 
while the bottom figure shows that the shear diverges less rapidly than the expansion scalar when approaching this initial 
singularity. As we discuss in section V, this presence of shear contributes to the enhancement of gravitational collapse and 
growth rate of large-scale structure in the Szekeres models. The amount of matter density (curvature) is varied as shown in 
the plots, and the spatially flat case is given for reference. RIGHT: Shear in negatively curved Szekeres Class-II models (or 
Class-I models with a fixed value of z). The same observations as for the positive case hold here. All four figures here also 
show that the amplitude of the shear increases as the amount of matter density does, so the shear is stronger in the positively 
curved models than in the flat case, which in turn is stronger than for the negatively curved cases. This is consistent with the 
respective strengths of the growth rate as plotted in figure [T] for the flat, positively, and negatively curved cases. 



non-linear terms into the differential equations governing the growth. This causes the enhancement of the growth of 
large-scale structure seen for the Szekeres models examined in the previous section. The growth was found not only 
to be stronger than that of the homogeneous Einstein-de Sitter model, but also stronger than that of the spherical 
(inhomogeneous) growth, making clear the contribution of the shear in the Szekeres models. 

Our plots for \a\ and as given by equations ([72)1 and ([7^)1 are shown in figure [3] for the flat and curved cases 
of the Szekeres Class-II models (the plots for Class-I with a fixed value of z are exactly the same, and so are the 
conclusions drawn) . The results show the time evolution well within the matter-dominated cosmological era and well 
after the radiation-dominated era. We chose a < 0.35, which is well before a w 0.60, the scale factor associated 
with equality between matter dominance and cosmological constant dominance in the FLRW-LCDM standard model 
[43|,[44|]. In other words, we are here only concerned with the cosmological era fully dominated by matter (dust). The 
upper plots show that the shear approaches infinity when a tends to zero (the initial singularity), while the lower plots 
show that the shear scalar diverges less rapidly than the expansion scalar when approaching this initial singularity, 
since the ratio is not diverging in that limit. These plotted behaviors for the exact scalars when approaching the 
initial singularity are consistent with the treatment of the asymptotic behaviors of the models derived by {26j using 
other considerations. The plots in the figure also show that the amplitude of the shear increases as the amount of 
matter density increases so that the shear is stronger in positively curved models than in the flat case, which in turn 
is stronger than the shear for the negatively curved Szekeres models. This is consistent with the respective strengths 
of the growth rate of large-scale structure as shown in figure Q] for the flat, positively curved, and negatively curved 
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FIG. 4: LEFT: Tidal gravitational field in positively curved Szekeres Class-II models (or Class-I models with a fixed value of z). 
The scalar \E\ and scalar ratio |-E|/p are plotted as a function of the scale factor well within the matter-dominated cosmological 
era (a < 0.35) (i.e., well before a « 0.60, the scale factor associated with matter and cosmological constant equality in the 
FLRW-LCDM standard model [43|, S3)- The plots show that the time evolution of the tidal gravitational field is indeed very 
similar to that of the shear, confirming the statements made in section V about the tidal field inducing shear in the matter 
flow. Similarly to the shear plots, the upper plots show that the tidal field approaches infinity when a tends to zero, while the 
lower plots show that the tidal field diverges less rapidly than the matter density when approaching this initial singularity. The 
amount of matter density (curvature) is varied as shown in the plots and the spatially flat case is given for reference. RIGHT: 
Tidal gravitational field in negatively curved Szekeres Class-II models. The same observations as for the positive case hold 
here. All the figures here also show that the amplitude of tidal gravitational field increases as the amount of matter density 
does, consistent with the amplitudes and evolution of the shear and the growth rate shown on previous plots. 



cases in the previous section. 

Again, it is worth noting that our first goal here is to represent the growth during the strictly matter-dominated era 
with a < 0.35. As shown in the plots, no curves experience divergences after the initial singularity, and all are far from 
any possible later time pancake singularity in the models f2fjj . This is similar to the plots of the energy density, which 
are certainly finite during the matter-dominated era of interest (see figure [2]). Additionally, such pan cake singularities 
can be avoided in some cases by the addition of a cosmological constant to this class of models [39[ . 

For further exploration, it is informative to consider the next propagation equation of interest (after the Raychaud- 
huri equation), and that is for the shear, given by (e.g., (2ll. l38l] ) : 

& a b + <J a c <J c b + ^6 a a b -±(5 a b + u a u b )a 2 = -E a b , (75) 

where E a t> is the mixed tensor corresponding the electric part of the Weyl conformal curvature tensor. The covariant 
tensor E ab is denned as (pll. l38j) 

E ab = C acbd u c u d , (76) 

where the Weyl tensor is given by 

& abed ^iabed + {gadRcb + 9bcRad ~ 9acRbd — gbdRac)/^ + R {9ac9bd — 9ad9bc )/6. (77) 

As discussed in [U [38|, E ab represents the tidal gravitational field, and its presence in the shear propagation 
equation (|75|) shows how the tidal field induces shearing in the source flow. The shear then affects the gravitational 
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collapse and clustering via the Raychaudhuri equation. Finally, it is well-known that in the Szckercs models the 
magne tic p art of the Weyl tensor, defined as H a j, = *C ac t,dU c u d , vanishes [THHH. (Here *C ac t,d is the dual of Weyl 
tensor (HIGH). 

Now, in order to analyze the contribution of the tidal gravitational field, we plot for the fiat and curved Szekeres 
Class-II models the time evolution of the invariant 



E = VE a b E b a . (78) 

For that, we first express the invariant in terms of the scale factor and measurable cosmological parameters as we did 
for the shear. We calculate the mixed components E a f, using equation ((75]) instead of the definition ((75)) . leading to 
the simpler expressions 

E 1 1 = E 2 2 = -2E 3 3 = + . (79) 

6 3 H 3 H a v ' 



It follows that 



2 / H H a \ 

E=± rAH +2 Ho • (80) 



Following the same steps as in the previous calculations for S, a, and 0, we find 



E = tV6^S. 81) 
a 6 



The definition of VLm allows us to write 



M_Om^ (g2) 



We can then get E in terms of G, H<j, and £l M using equations (|49|) and (|50[) . Equation (|81j) gives finally 



\E\ = ^-^E*G. (83) 

\E\ 



Similar to the shear over the expansion ratio, it is useful to define here the scalar ratio - — - (see, for example, [26j) in 
order to study early time divergences. Using the expression (|52|) for the energy density we can write this ratio as 

\E\ 1 aG 



p V6l + aG' 



(84) 



Our plots for the tidal gravitational field scalar \E\ and the ratio ^ are given in figure 2] for Class-II (the plots for 

Class-I with a fixed value of z are exactly the same, and so are the conclusions drawn). The scalar \E\ and scalar ratio 
\E\/p arc plotted as a function of the scale factor well within the matter-dominated cosmological era. We find that 
the time evolution of the tidal field is very similar to that of the shear and confirms the discussion above about the 
tidal field inducing shear (according to the shear evolution equation (|75p) in the source matter flow. We also see the 
same behavior as for the shear at earlier times where the tidal field diverges as a approaches the initial singularity, 

I -El 

but the ratio L-L does not, because the tidal field diverges less rapidly than the matter density as the scale factor 
goes to zero. The amount of matter density (curvature) is varied, and the figures show how the amplitude of tidal 
gravitational field increases as the amount of matter density does, consistent with what we observed for the shear and 
the growth rate plots. 



VI. CONCLUDING REMARKS 



We considered the growth rate of large-scale structure in the Szekeres inhomogeneous cosmological models written 
in the Goode and Wainwright representation. Using the Raychaudhuri equation, we derived exact equations for the 
growth rate for the two Szekeres classes. We explicitly expressed the equations in terms of the under/ovcrdensity 
and measurable cosmological parameters, leading to further insights on the growth rate of structures in the Szekeres 
models and also putting them in a framework close to comparison with cosmological observations. 



16 



In Class-II, the background density is only a function of time, so we defined the undcr/ovcrdensity in the standard 
way, while in Class-I we defined an invariant under /overdensity using the quasi- local averaged density over a spatial 
domain. For both classes, we found that writing these equations in terms of these under/overdensities instead of 
the metric functions allows for the growth equations to be remarkably split into two meaningful parts. The first is 
identical to the usual growth equations of a perturbed matter-dominated FLRW model. The second part is similar 
to second-order perturbations, but here the equations derived are all exact, and this part represents the nonlinearity 
of the Szekeres exact solution. 

We integrated numerically the exact equations obtained for flat Class-II, curved Class-II, and curved Class-I models. 
We note that flat Class-I models have no growing modes and so are of no interest to our purpose here. For curved 
Class-I models, we use a domain delimited by a fixed value of z. In all the cases, we found that the Szekeres growth 
rate is up to 3-5 times stronger than that of the perturbed FLRW models. The Szekeres growth (with a corresponding 
Einstcin-dc Sitter background) is also found to be stronger than the well-known nonlinear spherical collapse with the 
difference between the two increasing with time. This shows that the growth is stronger and different when we use 
the more general Szekeres inhomogeneous models where shear and tidal gravitational fields are present in the matter 
source. 

In order to explore this Szekeres strong growth further, we derived explicit expressions for the shear scalar and 
tidal gravitational field part of the Weyl tensor, again all in terms of the scale factor and measurable cosmological 
parameters. Our analysis shows and confirms how the shear acts in the Raychaudhuri equation like an "effective 
source" along with the energy density to enhance gravitational collapse and produce stronger growth of structure in 
the Szekeres models. We also derived and plotted the time evolution of tidal gravitational field which induces shearing 
in the cosmic fluid flow, as can be seen from the shear evolution equation. 

In this analysis, we focused our interest and results to be well within the matter-dominated cosmological era. 
That is well after the radiation-dominated era and also well before the equality time between matter dominance and 
cosmological constant dominance in the FLRW-LCDM standard model. In other words, we are here only concerned 
with the cosmic era fully dominated by matter. We also plotted the matter energy density during this era and found 
that it starts diverging (as expected) close the initial singularity, but after that it decreases monotonically with no 
other divergences during this interval, showing that we are far away from any possible later time pancake singularity 
in the models. Additionally, such pancake singularities can be avoided by the addition of a cosmological constant to 
the models as shown in other works. We also plotted the time evolution of other scalars, such as the shear over the 
expansion and the tidal gravitational field over the energy density, and found that they tend to zero when approaching 
the initial singularity, in agreement with previous results that used other considerations. 

It is worth noting that the enhancement of the growth found here in the Szekeres models during the matter- 
dominated era could suggest a substitute to the argument that dark matter is needed when using an FLRW-LCDM 
model during the matter-dominated era in order to explain the enhanced growth and all the large-scale structure 
that we observe today. Indeed, it is a well-known argument in FLRW models that dark matter is necessary to 
enhance the growth during the matter-dominated era in order to explain all the presently observed galaxy clusters 
and supcrclusters. In the Szekeres models the enhanced growth seems present with no requirement of dark matter. 
A similar conclusion was reached in, for example, [53| but from analyzing growth in a modified gravity model. Of 
course, inhomogeneous models with the presence of shear and a tidal gravitational field should also be explored to 
analyze galaxy rotation curves and other arguments in support of dark matter. We plan to explore this in follow-up 
work. 

Finally, the current results for the growth in Szekeres expressed within standard schemes of growth rate studies 
and in terms of observable cosmological parameters will be useful toward building a thorough framework where 
inhomogeneous Szekeres cosmological models can be compared to various current and future cosmological data sets. 
Comparison of the growth rate in Szekeres to available and future data probing the growth is out of the scope of the 
current paper but will need to be done in future and follow up works. 
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